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Abstract 

We present and analyze a multigrid algorithm for the acoustic single 
layer equation in two dimensions. The boundary element formulation of 
the equation is based on piecewise constant test functions and we make 
use of a weak inner product in the multigrid scheme as proposed in [4]. 
A full error analysis of the algorithm is presented. We also conduct a 
numerical study of the effect of the weak inner product on the oscillatory 
behavior of the eigenfunctions for the Laplace single layer operator. 

This paper is dedicated to Prof. G.C.Hsiao on the occasion of his 75th 
birthday. 

1 Introduction 

A model for the scattering of acoustic waves by a bounded obstacle is given 
by the Hclmholtz equation in the exterior of the scatterer, with appropriate 
growth conditions on the scattered field. One can reformulate this problem in 
terms of integral equations on the surface of the scattering object via direct 
or indirect boundary integral formulations. Or one may consider scattering 
from a coated bounded obstacle, in which case an integral equation can be 
used to prescribe a non-reflecting condition on an artificial surface surrounding 
the object. In both, one has to find numerical approximations to solutions of 
boundary integral equations. In this context, Galerkin type methods have been 
studied extensively and become popular over recent years, see for example the 
monographs [23], [16]. 

Our focus in this paper lies on integral equations of the first kind, which 
arise naturally in the direct boundary integral method for the Dirichlct problem. 
The main integral operator involved is called the single layer operator and may 
be viewed as a pseudo-differential operator of order minus one. Several authors 
have observed advantages of using integral equations of the first kind (e.g. [17]), 
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for example when the scattering object is very thin (e.g. [15]) or when the 
scattering surface is not smooth. Indeed, for problems of crack propagation in 
elasticity or scattering from a screen, integral equations of the first kind are the 
most appropriate model. 

Other popular integral equation strategics include the use of combinations 
of single and double layer operators, to avoid issues of non-uniqucness and the 
potential for numerical instability near possible eigenvalues of the operators. 
These approaches include the famous Brakhage- Werner and Burton-Miller for- 
mulations. In these cases as well, discrete strategies rely on effective and accu- 
rate approximation methods for the layer operators involved. 

Due to the non-local behavior of boundary integral operators, they typically 
lead to dense linear systems upon discretization. Though one only needs to 
mesh on a surface of co-dimension one, the fill-in in the matrices corresponding 
to the integral operators is significant. Without some form of preconditioning 
or acceleration, these methods then become prohibitive. 

One possible preconditioning strategy is the use of a multigrid scheme. How- 
ever, the use of standard multigrid smoothing operations is inappropriate for 
negative order pscudodiffcrcntial operators. Such operators link highly oscil- 
latory eigenfunctions to small magnitude eigenvalues. This ruins the basic in- 
terplay between standard smoothing of oscillatory error components and the 
possibility to represent the remaining error components on coarser grids. The 
remedy for this is the use of a weaker inner product in order to modify this 
spectral feature of the operator. This approach has first been described in [4] 
for positive-definite operators. We follow the same path for the acoustic single 
layer equation. The numerical examples in Section 4 exemplify how the spec- 
tral behavior of the discretized operators changes through the use of the weaker 
inner products. 

Several related works can be found in the literature. The single layer equa- 
tion associated to Laplace's equation has been treated and analyzed in [12] using 
a BPX preconditioner. The same equation was considered in [19]. Here, the 
authors studied a multigrid method for large-scale data-sparse approximations 
to the single layer operator. In the acoustics case, the use of Haar basis func- 
tions and compression type multilevel algorithms for integral equations of the 
first kind has been studied in [21]. Algebraic multigrid preconditioners, based 
on the smoother in [4], have been developed in [18]. 

The purpose of this paper is to prove convergence of the multigrid algorithm 
given in [4] when applied to the indefinite acoustic single layer discretization. 
We will make use of perturbation- type arguments such as in [3] and [14]. We 
also state the algorithm in the situation where a non-uniform discretization is 
used. 

The design of the algorithm heavily relies on the above references and we 
have reported its promising numerical performance in [13]. Currently, our codes 
assemble and multiply matrices in 0(n 2 ) complexity, and so the computational 
cost of the iterative solution by multigrid is the same. The efficient assembly and 
matrix multiplication of the corresponding matrices is an active research issue in 
multipole and hierarchical matrix theories. Our numerical results indicate that 
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we can approach optimality in complexity through multigrid, once assembly and 
matrix multiplications are done optimally. 

At this juncture, we note important directions for extending this work, moti- 
vated by recent developments in the study of boundary integral equation meth- 
ods. Firstly, the use of combined integral equation formulations, including those 
by Brakhagc- Werner [1] and Burton-Miller, provide stable formulations for the 
solution of scattering problems. These are popular approaches, and the devel- 
opment of a multigrid strategy for these would be of interest. Both the methods 
of implementation and analysis would be different than those in this paper, to 
account for the different properties of the combined layer operators. Secondly, 
the methods presented in this paper are not tailored to high-frequency scat- 
tering, and our analysis does not include wave-number explicit bounds. Wave- 
number explicit bounds are described, for example, for combined field operator 
approaches [8, 9, 10, 22, 20]. The scale resolution condition of kh being suffi- 
ciently small (for piecewise linears) needs to be met; in high-frequency settings 
non-polynomial approximation spaces may be a better approach, [7]. The de- 
pendence of the coarsest mesh on the wave number is assumed to satisfy this 
requirement. In this paper, our focus is on a much simpler situation: how to 
use a first-kind integral equation with piecewise constant approximants to the 
solution, in the low to medium frequency situations of scattering from polygonal 
domains. 

We now give a brief derivation of the acoustic single layer equation using the 
framework of a direct boundary integral approach. We consider the following 
exterior Hclmholtz problem with prescribed Dirichlet data on the boundary of 
a scatterer. Here, T is a simple, closed polygonal curve in the plane and VL ext 
denotes its exterior domain. 

—Ait — k 2 u = in Q ext , u = g on T and lim (— inu) = 0. 

r— > oo Or 

To guarantee unique solvability we assume a non-zero, real wave-number n G 
K, such that k 2 is not an interior eigenvalue of —A. The Sommerfeld radiation 
condition is given in terms of the usual radial component r in polar coordinates. 
It is well known that the solution to this problem is fully determined by its 
complete Cauchy data g — -f + u and a = B+u, where j + : Hl oc (U, ext ) — > 
H l / 2 {£) and B+ : H^ oc (n ext ) -> H~ l / 2 (T) denote the exterior trace operator 
and the exterior conormal derivative respectively. The normal n is assumed to 
be outward to $l ext , ie, it points into the bounded region. In fact, for x € Q ext 
one has an integral representation formula for the solution to the boundary value 
problem (e.g. [16], [24]), namely 

<*) = -{ j r 4 1 \n\ X -y\)a(y)d Sy + j jf = Vl) g(y)ds y . (1) 

The kernels of the two integrals are given in terms of the Hankel function (z) 
and its conormal derivative. According to the representation (1), it is sufficient 
to find the unkown surface density a . To do this, one exploits the jump relations 
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of the two integrals, i.e. their behavior in the limit as x approaches V from both 
the interior and the exterior of the scattering domain. These relations appear 
when we take the trace of equation (1) and lead to the following integral equation 
for en 

Va = f G H 1 ' 2 {Y). (2) 

Here, the right hand side / = (\l + K) g depends on the Dirichlet trace g 
and requires the evaluation of the double layer operator K. The single layer 
operator V and the double layer operator K are both defined in terms of singular 
integrals. 

V : i?- 1/2 (r) -> H 1 ' 2 ^), Va(x) 

K : H X / 2 {Y) -> H l ' 2 {Y), Kfi{x) 

One should note that the above approach to reformulate the original bound- 
ary value problem into an integral equation on Y is by no means unique. How- 
ever, several methods lead to a single layer equation of the form (2). In the fol- 
lowing, we are interested in the weak form of equation (2). Given / G i? 1 ' 2 (r), 
find a G H~ l / 2 {Y) such that 

V(a, n) = (/, fi) for all M G H-^ 2 (Y), (3) 

where the continuous sesquilinear form V : H^ 1 / 2 (Y) x iJ -1 / 2 (r) — >• C is defined 
by V(a, fi) = (Va, fx), and (•, •) denotes the duality pairing between if 1 / 2 (r) and 
i? _1 / 2 (r). The single layer operator which corresponds to the Laplacian will be 
denoted by A. 

A : H-^ 2 (T) — ► H^ 2 (T) 7 Aa(x) := / ln(\x - y\) a{y)ds y x G Y. 

We note that the underlying differential operator is the principal part of the 
Hclmholtz operator. Its associated sesquilinear form A(-, •), denned similarly to 
V(-, •) above, is positive definite after the region has been scaled properly (see 
[24] for details), i.e. there holds 

A(a,a)>C\\a\\ 2 H . 1/2(r) for all a G H^^iY). (4) 

Consequently, it defines an inner product whose induced norm is equivalent to 
the natural energy norm in iJ _1 / 2 (r). This will play an important role in the 
analysis of Section 3. 

The paper is organized as follows: In Section 2 we present a multigrid al- 
gorithm for integral equations of the first kind, and introduce a (computable) 
inner product whose induced norm is equivalent to the natural norm in H^ 1 (Y) 
on finite dimensional test spaces. The multigrid strategy relies on reformulating 



:= l -jH^{K\x-y\)a{y)ds y x G Y, 

i f dH { 1} ( K \x-y\) 
:= 4j r ^ Ky)dSy XeT - 
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Problem (3) using both the standard inner product in £r -1 (r) and the new 
computable version for piecewisc constant functions. We introduce smoothers, 
and present a matrix version of the algorithm. Section 3 consists of a conver- 
gence analysis. The key component is a careful study of the difference between 
the single layer operators for the Laplacian and for the Helmholtz equation. 
We conclude this section with a convergence result. Finally, in Section 4 we 
present some numerical experiments describing the spectral behaviour of the 
single layer operator, as well as that of the operators used in the reformulated 
problem. We see, in the context of a smooth curve (a circle) and a Lipschitz 
curve (a square) how the use of the weaker inner product renders the problem 
suitable for a multigrid strategy. We conclude by reporting on the convergence 
behavior of the algorithm applied to simple test cases. 

2 A multigrid algorithm 

In this section we present the multigrid algorithm, originally proposed in [4] for 
the positive definite pscudodifferential operators of order minus one, and which 
was applied to the acoustic case in [13]. The multigrid algorithm is presented 
below in terms of a smoother, whose definition is postponed to Subsection 2.2. 
The smoother is realized using a weak base inner product. We make this more 
precise in Subsection 2.1. 

Wc first establish notations needed to describe the multigrid algorithm. As- 
sume that the polygonal boundary T is composed of finitely many straight edges 
Tj. Each Tj is meshed by a coarse grid of line segments of length Zj. 

We successively refine this grid in a uniform way by breaking each element 
in half and adding the respective midpoints to the vertices of the next finer 
level. On every level of refinement k — 1, 2, . . . , J, we label the vertices in such a 
way that x\,x\,..., x k Nk , x k Nk+1 — x\ is a counterclockwise enumeration. Now, 
let cj) k be the characteristic function of the line segment if — conv(x^, 
(i = 1, 2, . . . , JVfc) and denote their span by Mk ■= span {</>f }■ For the sake of 
easier notation we will suppress the level number k in our notations whenever 
the context rules out any ambiguity. This construction yields a sequence of 
nested finite-dimensional spaces 

AiicM 2 c...cMjC H~ l ^(T). 

We now define the discrete operators Vk ' Mk — ► Mk with the help of the 
i? _1 (r) inner product, denoted by (■, The defining relation is 

(Vfecr, = V(er, fi) for all a, fi £ M k - (5) 

Analogously, we choose fk € M k to satisfy (fk,n)-i = (/, /i) for all [i G M.k- 
Then, on every level k, the equation of interest can be written in operator form 

as 

V k cFk = .fk- (6) 
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In order to describe the algorithm in a function space setting, we shall also need 
the H^ 1 - projections Q k ■ iJ _1 (r) — > A4 k , which are defined by 

(Qfccr, /u)_i = (cr, forall^eA^fe. 

We further need a family of smoothing operators R k ■ Mk — > Mk- It is defined 
precisely in Subsection 2.2, but for now, we just assume that R k are some given 
linear operators. Then, given an initial guess o~q € Mj, the multigrid iteration 
computes a sequence of approximate solutions to (6) using an iteration of the 
form Oi+i = Mgj(ai,fj), where Mgj(- , •) is a mapping of Mj x M. j into M. j, 
defined recursively by the following algorithm: 

Algorithm 1. Set Mg 1 (a, f) = V\~ x f . If k > 1 we define Mg k (<r, f) recursively 
as follows: 

ax=a + R k (f-V k a), (7) 
Mg k (a,f) = a 1 + Mg^i^Qk-iif - V k a x )). (8) 

This is a simple variant of a V-cycle multigrid scheme, which only uses pre- 
smoothing. Equivalcntly, we can write the iterative scheme as a linear iteration 
method 

o-i+i = o~i + Bj (/,/ - Vj (Tj), 
with an "approximate inverse" Bj : Aij > Aij defined by 

B k fk = Mg k (0, f k ) for all f k e M k and k = 2, . . . , J. 

This operator is useful as a preconditioner in preconditioned iterative methods. 
The matrix version of Algorithm 1 is given in Subsection 2.3. 

2.1 Discrete inner products 

The use of the i/ -1 (T) inner product in defining the operators V k and Q k for 
the multigrid algorithm confronts us with the question of computability. We 
will have to work around this issue by introducing equivalent computable inner 
products. These inner products will be used to define smoothers in Section 2.2. 

The problem of calculating the i/ _1 (r) inner product of two elements in M. k 
is related to the solution operator of a second order boundary value problem on 
the boundary curve, namely 

-u" + u = v. (9) 

Here, v G i? _1 (r) is a given function, u has periodic boundary conditions, and 
the primes denote differentiation with respect to arc-length. The weak form of 
this problem is uniquely solvable and we denote the bounded solution operator 
by T : H' 1 ^) — > H 1 ^). Then, for v,w € H' 1 ^) it is easily verified that 
(v,w)_i = {Tv, w)t = (v,Tu>)r, where (■, -)r denotes the complex L 2 (r)-inner 
product. Unfortunately, the use of the exact solution operator T is infeasible. 
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Instead, we discretize (9) using a second-order finite difference method. This 
finite difference method results in an N k x N k linear system with an tridiagonal- 
like matrix k k . We need to introduce some more notation in order to see the 
details. Functions in Mk can be represented through their basis expansions 
with respect to the (f>^. This is done via the following map: 

e k : M k — > £ N \ [e k (<j)] t = 1 (a,^ fc ) r , (10) 

H 

where Zf := meas(r 4 fc ). Since the basis functions <\>\ are the (orthogonal) in- 
dicator functions of the segments r*, the basis expansion for any c 6 A^t is 
a = Yli=i [ e k(c)]i ) so the map in (10) gives the vector of coefficients. 
We then define the (invertible) operator A k ■ Mk — > Mk through 

A „ „ ST f i e ( a )h+ - i e ( a )h [e(cr)] t - [e(cr)] t _ ^ fe 

Aka = a q hi~ ) 

where we use the notation i+ = (i + 1 mod Nk) and i— = (i — 1 mod N k ). 
As in (11), we will often drop the superscript k and identify l\ and k to be the 
same for convenience. 

The inverse operator A^ 1 or equivalently the inverse matrix A^ 1 serve as 
approximations for the solution operator T. This motivates the definition of 
the computable, discrete inner products on the spaces Mk via 

[faiflk := (AfcV^) r f or all 0,^e M k . (12) 

The following lemma shows that || • ||_i and [•, are equivalent norms, with 
the equivalence constants independent of the refinement levels k = 1, J. This 
result can be found in [4] and [12] for the case of the screen problem when V 
is in fact an open boundary patch or an open line segment respectively. Here, 
we give a proof for a closed boundary curve T. For the analysis we will assume 
that all the mesh element lengths if are such that C\h k < l\ < Cih k for 
some fixed positive constants Ci, C2. Here hk is a representative mesh size, e.g., 
hk = max(meas(Tj fe )). 

For the proof we need a standard approximation result. Let 9 £ H 1 ^) and 
let 9 k € M k denote the piecewise constant approximation defined by 

9 k :=Y^6(x*)4>l (13) 

This is well defined for 9 in H 1 ^) by the Sobolev inequality 

Ph.^r) < C \\9\\ HHT) for all 9 G H\Y). (14) 

Then 

\\G-9 k \\ LH r)<Chk\9\ H i (T) for all 6 e H X {Y). (15) 
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To prove this, it suffices to observe that 



dx 



\\e-o k \\i. 2[Tk) = / \9(x)-e(xt)\'dx-- 

< j Tk (/r i|0,(e)|2 ^) ( x ~ x i) dx 

= f (/,iw«)=|i*iWw 

Equation (15) follows by summing over all elements. We will use (15) in the 
proof of the next lemma. 

Lemma 1. There exist constants c, C > independent of J such that 

CM 2 .! < [<7, ( r] fe <C'|| ( r|| 2 _ 1 /or a g M k (fc = l,...,J). (16) 
Proof. Given an element 

N k N k 

a = J2[e{a)]i we define a = ^I e ( cr )]i 

where V'f is the continuous and piecewise linear function which takes the value 
one at node xf and vanishes at every other node (otherwise known as the "hat" 
function). Note that a is in H 1 ^). Then applying (15) to a, we obtain 

Ik - ^!U 2 (r) < Ch k ||<9<7|| L 2 (r) , 

where da denotes the derivative of <r with respect to arc length. Note that 

^ = E Ae(^ + -[e(a)]A , ^ 

i ^ ' 

From this it follows by straight forward calculations that 

\t\m { r)<Ch^ 2 \\o-\\l Hr) . (18) 

Before proving the inequalities, we make a few observations. Using the 
operator from (11) we see for a, /i <E M. k , 

(A k a, fj,) r 

= i^Nr -i^J y j2 & & J i e WU<i>i ds 

= Mr - £ ~ [e(g)] * ~ [e(g)] \ [e(g)] " ) Mi 



2 



(a, /i)r + / dadfsds. 
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by (17). The right hand side above defines a sesquilinear form 

a(a,p) = (<7,n) r + (<9ct, <9/2) r . 

It is easy to see that 1 1 c 1 1 C r ) ~ lkl!i 2 (r)' w here X ~ Y indicates that C\X < 
Y < C2Y holds with some positive constants C\ and C2 independent of the 
mesh. Thus we have proven that 

||5ili2-i(n ~ a(cr, a) = {A k a, a), for all a G Mk- 

If we combine the above norm equivalence with (18) and write ||o"||^ = (A k cr, a), 
it immediately follows that 

h\\l k =a(a,a) = \\a\\l 2{r) + \a\ 2 HHr) < (1 + C h^) ||a||| 2(r) < C hf \\a\\ | 2(r) . 

Here, C denotes a generic constant independent of mcshsize. In other words, 
this yields the inequality 

A max (A fc ) < Ch^ 2 , (19) 

which could also be shown using other methods. 

Wc now begin proving the inequalities of the lemma, starting with the second 
inequality in (16). 



(A k a,a) r = a(a,a) = sup — ' < C sup M _ ll2 

flEM k M6M t Gi||/z|| ffl(r) 

, n |(A fc a,/2)r| 2 + |(^ fc a,M-/2)r| 2 

i o sup 1| — lis 

fieMfc ll/ x llj?i(r) 

|(^fc(j,^)r| 2 ll^^ll|2 ( r)ll^-All|2(r) 

< G SUp — 1| _„2 1 11 ~||2 

neM k llMlliji(r) HWIif^r) 

„ |(AfcCT,^)r| 2 „ ^ 2 n . || 2 

< G sup 1 p h G fr fc ||A fe <T|| L2(r) 



= C (\\A k a\\ 2 H - 1{r} + h 2 k \\A k a\\i Hr) ) . (20) 
The inverse inequality 

llMlU^r) < Chi 1 \\n\\ H -i(r) for all A* G -Mfc (21) 
can be found in [4, Eq. (3.20)]. Using it in (20), we get the upper bound 

{A k a,a)v <C\\A k a\\%- 1(r) for all a G 7W fe (22) 

or equivalently 

(AfcV.M)r<C||A*||| r - 1( r) for all /i G .M fe . (23) 
This is the upper inequality stated in the lemma. 
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It remains to prove the first inequality in (16). For this, we need a stability 
result for the i 2 -orthogonal projection n k '■ H 1 ^) — > JA k , namely 

(A k n k e,n k 6)<c\\e\\ 2 m(r) . (24) 

The result (24) follows once we prove that for all 9 G i/ 1 (r), there exists 9 k € 
Mk, such that 

\\0k\\A k <C\\9\\ HHr) (25) 
\\0-0k\\mr)<Chk\e\ HHr) . (26) 

Indeed, writing II k 9 = II k {6 - 6 k ) + 6 k , 

\\n k e\\ Ak < \\n k (6-e k )\\ Ak + \\e k \\ Ak 

< ch- 1 \\n k (e - e k )\\ L 2 {r) + \\e k \\ Ak by (19) 

<Ch^\\(e-9 k )\\ L 2 {r) +C\\e\\ m(r) by (25) 

<C||0||jji ( r) by (26). 

Therefore, let us now exhibit a 9 k satisfying (25) and (26). 

Consider the 9k defined in (13). By (15), we get (26). That 9k also satis- 
fies (25) is seen using the L 2 -orthogonal projection Q k into the space of continous 
functions on T that are linear on each r*. It is well known that Q k is stable in 
the i? 1 (r)-norm [6]. Hence, using the standard approximation properties of the 
projection Qk9 and the linear interpolant 9 k of 9, we find 

\0 k \ H \r) < \k ~ Qk9\ H ^v) + \QkO\ H ^v) 

< Ch^ \\9 k - Q k 9\\mv) + \QkO\ H i(T) 

< Ch^(\\9k - 6\\ LHT) + \\9- Q k 9\\ L 2 {T) ) + \Qk9\m { v) 
<C\\9\\^ {T) . (27) 

Moreover, by the Sobolev inequality (14), we also have 



INI 2 = £|0(*i)| a IWI a ^ IMIW) (E/ X ) < C Pfm(vy 

Therefore, combining (27) and (28), we have 

\\9k\\ Ak = ll^fclL=(r) + \Qk\m{T) < C ll^lliji(r)> 

which proves (25). 

To complete the proof of the first inequality in (16) of the lemma, 

a M) 2 (<J,n k 9) 2 
wWh-hf) = SU P jTTTp = SU P T^p 



(28) 



eefl-i(r) \A- k ii k v,iikV) 



where we have used (24). 
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2.2 Smoothers 

Using the inner products from Section 2.1 we can define a simple Richardson 
smoother suitable for multigrid algorithms (see for example [2]). The smoother 
is given by 

[ii fc M]fc = J-M)_x, (29) 
Afe 

where X k is related to the Rayleigh-Ritz quotient involving the positive definite 
sesquilinear form A(- , •), namely 

A, = sup f^l, (30) 

as follows: We assume for our analysis that Xk is a number that satisfies 

X k <X k < CX k (31) 

for some mesh independent constant C. 

Define the operator A k as in (5) but with the sesquilinear form V(-, •) re- 
placed by A(-, •). The eigenvalue X k is then a computable version of the largest 
eigenvalue of the operator A*, with respect to the minus one inner product. In 
practice, we could choose X k = Xk, or an approximation to X k computed by a 
few iterations of the power method. 



2.3 Matrix version 

Now we give a readily implcmentable matrix version of the previously given 
multigrid algorithm. Since M. k C AAk+i we can find numbers c^; such that 
$ = Y^i 1 ( t>'i +1 ■ These entries define the N k x N k +i restriction matrix 
Cfc by [Cfejj,; = Ci i. This matrix and its transpose are used as intergrid transfer 
operators. We further define the operator 

f k :M k ^C Nk , [fk(o-)) i = (o-,$)- 1 . (32) 

Algorithm 1 can then be translated into an approximation scheme for the matrix 
version Vju = b of equation (6). Here, the vectors are given by b = f k (fj) 
and u = e k (a,j) and the system matrix is Vj = [{Vj(f>j , <j>(})i,j- This yields 
a procedure Mg 7 (s,b) that outputs an approximation to the solution given an 
input iterate s. To describe it we will also need a matrix of the operator A k , 
which we denote by Afc. Specifically is the matrix satisfying e k (A k <r) — 
Afc e k (o~). It is a circulant matrix with cyclically shifted rows of the form 

A fc = Circulant [0 • • • 0, -(Mi-) -1 , 1 + if + {h k-)' 1 , -if, • • • 0] . 

Note that A& is neither tridiagonal nor symmetric, but E k A k is symmetric, where 
Hfc is a diagonal matrix whose i th diagonal entry is meas^). The translation 
of Algorithm 1 into its matrix version is done via the identities of the following 
lemma. 
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Lemma 2. The following identities hold: 



fk (V k g) 


= V fc e fe ( ff ), 


for 


all g e Mk, 


(33) 


ffe-i(<9fc-isO 


= C fe _i ± k (g) 


for 


all g £ Mk, 


(34) 




= C t fe _ 1 e fc _i( ff ) 


for 


all g £ 


(35) 


e k (R k g) 




for 


a// ,g g A^fc, 


(36) 


eiiV^g) 




for 


all g e Mi. 


(37) 



Proof. Let us prove (33): 

ftC%<M = (V fc ff,#)-i = V(ff,#) 

N k N k 

= X)[e*(fl)]iV(^,^) = E[ V *k,Mff)b 

3=1 3=1 

= [V fc e fc (.g)]j. 

Next, let us prove (34): 

N k 

%-x{Qk-xg)]i = (Qfc-iS,^ 1 )-! = (<?,^ _1 )-i = (ff,X>k-iki#)-i 

z=l 

= [C fc _i4(fl)]i, 

since Cfe_i is real. The proof of (35) is similar. To prove (36), observe that 

i[J*(ff)]i = i(ff,0?)-i = [i?fcS,^]fc = M^efc^s)]*. 

Multiplying both sides by the symmetric matrix A^H^ 1 , we obtain (36). Proofs 
of the other identities are similar. 

These identities enable us to state a matrix version of Algorithm 1. For 
example, applying ej, to the step (7) of Algorithm 1 and using (36) and using 
Lemma 2, we have 

e fc (tri) = e k (<j) + e k {Rk{f -VkO-)) = e k {<j) + \- k 1 K k - 1 kif k {f -V k <j)) 
= e k (a) + A^ 1 H^ 1 A* (f fc (/) - V fc e fc (tr)). 

Thus, the matrix version of this step is si = s + A^" 1 1 k k (b — V^s) with s\ = 
e/ c (ai),s = efe(cr), and b = ffe(/). Using also the other identities in (33)-(37), 
we can similarly translate the entire algorithm. We then obtain the following 
matrix version of the algorithm Mgj(s,b), which outputs an approximation for 
the solution of the matrix equation V,/u = b, given an input iterate s. 

Algorithm 2. Let s and b be any given vectors in C Nk . Define Mg fe (s,b) 
recursively as follows. Set Mg 1 (s,b) = Vj^b. If k > 1, define Mg fe (u,b) as the 
vector in C Nk obtained recursively by: 

si = s + A,: 1 H,: 1 A^(b-v fc s), 

Mg fe (s, b) = si + Cj^Mg^O, C fc _i(b - V fc si)). 
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It is important to note that the inverse of is not needed in the implemen- 
tation. We only need to multiply by k\. In the case of a uniform mesh with 
mesh size hk the multiplication by the matrix H^ 1 reduces to multiplication with 
the constant 1/hk- Furthermore, it has been shown in [4] that 

X k = 0(l/h k ). (38) 

Motivated by this observation, we could replace the expression X^ 1 1 in Al- 
gorithm 2 by C, for some constant C (which depends on the geometry of the 
domain). This would allow us to bypass the eigenvalue computation, which 
would otherwise be inherent in the algorithm. The numerical experiments we 
present later, however, explicitly include the expression A^T 1 . 

Note that a matrix preconditioncr for is implicit in Algorithm 2 and is 
defined by B^b = Mg fc (0,b). A theoretical study of the convergence rate of the 
algorithm is presented in the next section. 

3 Convergence Analysis 
3.1 Preliminary Steps 

Before we can give a detailed description of the convergence behavior of Al- 
gorithm 1, we need to pave the way with some preliminary discussions. In 
particular, we need the Galerkin projections Pk : i/~ 1 / 2 (r) — >■ Mk satisfying 

V(P k a, (j,) = V(cr, fi) for all /i G M k - (39) 

As we see next, such operators are well defined once the mesh size is sufficiently 
small. The assumption on the wave number that n 2 is not an interior eigenvalue 
of —A implies that for homogeneous right hand side the equation Va = 
only has the trivial solution a = 0. It is then a standard theorem on compact 
perturbations (see for example [23, Theorem 4.2.9]), that the discrctized version 
of equation (2) has a unique solution a k G M.k if the corresponding meshsize 
hk is sufficiently small. Furthermore, wc know that the Galerkin solutions ak 
converge quasi-optimally to the true solution <r, i.e. 

||er- Cfc || ir-i/2(r) < C min II "- «fc|U-i/afrv 

Also, once the mesh size is sufficiently small, the Galerkin solutions depend 
continuously on the data, i.e., 

K||ff-l/2 (r) < C \\f\\ H1 ,2 {T) . (40) 

As a consequence, wc immediately have the following lemma which shows that 
Pi is a well defined continuous operator once hi is small enough (and so is Pk 
for k > 1). In the lemma and elsewhere, wc write || • ||a for the vector norm 

1/2 

A(- , •) . We will also use the same notation for the operator norm induced by 
this vector norm. 
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Lemma 3. There exists an H > such that once the coarse mesh size h\ is 
less than H , there is a unique PkO~ satisfying (39) for all k > 1 and moreover, 

II^HU < C ||<t||a. (41) 

Now let us introduce a few ingredients needed to analyze the multigrid al- 
gorithm. A simple induction argument shows that Mg k (- , •) as defined in Algo- 
rithm 1 is linear as a mapping from Aik X M-k into AA k - It is also consistent in 
the sense that Ok = Mg k (ak, VkCTk) holds for all o~k € .Mfc. The error reduction 
operator of the scheme is given by 

E = M flj (.,0) > (42) 

i.e., if e* denotes the error at step i, we have e 8+1 = Mg J (e l , 0). Furthermore, the 
error reduction operator admits a product representation as shown in Lemma 4. 
This representation will be essential in the convergence analysis of the V-cycle 
scheme. Proofs of such results can be found in [2]. 

Lemma 4. Let Tk = RkVkPk for k > 2 and set T\ = P\. For k > 1 we then 
define EkU = u — Mg k (0, VkPku) and set Eq = I , the identity operator. Then, 

E k = E k _ 1 (I-T k ), and (43) 
E = (J — Ti)(J — T 2 ) •••(/ — T 7 ). (44) 

The same representation holds for the error reduction operator E of the def- 
inite problem. Analogous to (39), we can define P k as the orthogonal projection 
into Mk with respect to the A(-, -)-inncr product. This is the Galerkin projec- 
tion for the principal part of the differential operator. If we set Tk = Rk^ k Pk, 
we get 

E = (I-f 1 )---(I-f J ). (45) 

This operator is proved to be a reducer in [4] . Specifically, in [4] the convergence 
for the symmetric version of the multigrid algorithm applied to the positive 
definite problem was shown. In fact, it was shown that the symmetric error 
reduction operator E s in this case is bounded away from 1 independently of the 
number of levels of refinement. The symmetric version differs from Algorithm 1 
by an additional post-smoothing step. However, it is well known (see, e.g., [5, 
Remark 3.4]) that the analogous result holds for Algorithm 1 with just the 
pre-smoothing, i.e., we have the following theorem. 

Theorem 1. The error reduction operator E for Algorithm 1 applied to the 
positive definite problem satisfies 

||E|| A <5<1, (46) 

where 5 is independent of J . 



14 



In order to analyze the algorithm for the indefinite Helmholtz case we look 
at the difference between E and E. Let = 1\. — Tk, and suppose for some 
positive a we have 

WZkU^dK for k = 1,..., J. (47) 

With this assumption, by well known arguments in an abstract multigrid set- 
ting [2, 14], we have the following theorem. 

Theorem 2. Let E satisfy (44) and E satisfy (45). Assume that (47) holds. 
Then, there exists a positive constant C 2 depending on C\, hi, and a above, 
such that: 

||E|| A < \\E\\ A + C 2 h?. (48) 

Wc know that ||E||a < 5 < 1 by Theorem 1. Hence by virtue of Theorem 2, 
to prove a convergence result for our multigrid application, we only need to verify 
the hypotheses of Theorem 2, namely (47). This will be done in Subsection 3.2. 

Before concluding this subsection, we need to establish one more ingredient 
for the multigrid perturbation argument. It is well known that the difference 
between the single layer potentials of the Helmholtz and the Laplace equations, 
namely D := V— A, is compact as a map i? _1 / 2 (r) — > iJ 1 / 2 (r). For our purposes 
we need the following lemma. 

Lemma 5. D is bounded as a map i/~ 1 / 2 (r) — > H 1 ^). 

Proof. In this proof, we use the explicit integral representation of the single 
layer potentials as given in Section 1. The operator D = V — A generates the 
sesquilinear form 

D(p,a) = (Dp, a), 
and is an integral operator whose kernel consists of the function 

f(x,y)=g{\x-y\) where g(z) = 1 -H ( ^\k,z) + — ln(z). 

The function g has the following asymptotic behavior as z approaches 0: 

g(z) ~ Cl + 0(z 2 log z) (49) 
g'(z)~c 2 (z\ogz) + 0(z) (50) 
5 "(z)~c 3 + 0(logz), (51) 

for some constants ci (depending on k). 

Let us now estimate the H 1 -norm of Da. Denote by df the derivative of / 
with respect to arc length along T. Then 

\\Da\\ 2 H1{T) = ||£>a||| 2(r) + \\d(Do)\\l 2{T) . (52) 
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Let us start by bounding the first term on the right hand side. Letting F x (y) 
f(x,y), we have 



\Do\ 



By the trace theorem 



L 2 (r 



f(x,y)a(y) ds y 
— J ll-^x||fl-i/2(r)||c||^_i/2( r ) ds x 



(53) 



for some ball B r of sufficiently large radius r (so that B r contains T) . The term 
1 1 -fa; 1 1 /^(B,.) is finite because 



^ y F x = g'(\x - y\) 



x-y 
\x-yV 



(54) 



is a square integrablc function (due to (50) and the boundedness of (x — y)/\x — 

y\)- 

By a change of variable (mapping x to 0), integrals of F x can be converted 
to integrals of Fq on transformed domains. Hence, by enlarging the transformed 
integration region, we have 

11^11^(5.) < \\F \\m {B2r} . 

This shows that the first factor in the integrand of (53) admits a bound inde- 
pendent of the integration variable x, so 



||-Do-||i,2(r) < j ll^llHi/2 ( r)!l cr llff-i/2(r) ds ^ 

< meas(r) C \\F \\ H i iB2r) \\<7\\ 2 H _ 1/2{r) 

ff-!/2(r)- 



<C||a|| 2 



(55) 



We treat the second term in (52) similarly. Denote by t x the unit tangential 
vector to T in the point itT, which is defined everywhere except on a set of 
measure zero (the corners). Then, 



\8Da\ 2 



where 



1= Jr 1 ^ X,V ^ a ^ dSy ' ' tx 

(^xf(x,y)) -t x a{y) ds y 

ffi/2(r) II" llH-V2(r) ua x 
G x {y) = t x ■ V x f{x, y) = -t x ■ V V F X . 



ds x 
ds x 



— I I II if i/2 epi II fll t/-i ii it-\ ds 



(56) 
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Note that by differentiating (54), 



V,G. = -g"{\x - y\) {x V) \l* {x — y) — g'{\x - y\)V v ( {x v) '** \ . 

\x-vr \ \ x -y\ J 

The term V y (x — y) ■ t x /\x — y\ is 0(l/\x — y\), while the multiplying factor 
g'(\x — y\) is 0(\x — y\ log \x — y\) by (50). Hence the last term is 0(log \x — y\). 
The first term on the right hand side is also 0(log |a; — y\), because of (51). 
Consequently, \7 y G x is locally square integrable on R 2 . Therefore, returning 
to (56), we can complete the estimation using a trace inequality and bounding 
1 1 G:r! I #1/2^) independently of x as before. Thus 

\\dDa\\l 2{r) <C\\a\\ 2 H . 1/Hry (57) 

Using (57) and 55 in (52), the proof is finished. 



3.2 Convergence 

Now we give our main result on the convergence of the multigrid algorithm for 
our application. The proof proceeds by verifying the hypotheses of Theorem 2. 
For this, we need a regularity result. Consider the solution s of the adjoint 
problem 

V(v, e) = F{rj) for all tj e iJ^ 1/2 (r), (58) 

for some linear functional F on i? -1 / 2 (r), or in other words F is in i/ 1 / 2 (r). If 
F is more regular, then we expect the solution e to be more regular. 
To make this precise, note that (58) can be rewritten as 

V(e,v) = F( V ), 

where V(-, •) is defined for smooth tr, fi by 

V(ct, m) = j^J^ ^H^\k\x - y\) a(y) n(x) ds y ds x . 

This form extends continuously to i7 _1 / 2 (r) x i/ _1 / 2 (r) and the operator V* : 
iJ" 1 / 2 (r) ^ i7 1 / 2 (r) defined by (V*a,fi) = V(a,/j,) is continuous. It can be 
written as 

V* = A + D* 

where D* is an integral operator analogously to D, but with an integral kernel 
conjugate to that of D. The same type of arguments as in Lemma 5 show that 

D* : if _1 / 2 (r) i y ff x (r) 

is continuous. Now, it is well known [23, Thm. 3.2.2] that for the positive 
definite problem Au = F, there is a regularity result: 

N|jja ( r) < C\\F\\ Hs +i {r) for < s < s 
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where sq is a positive number depending only on the angles of T. Applying this 
result with s = to (58) rewritten as V*e = F, or in other words, Ae = F — D*e, 
we obtain that 

\\s\\ L . {T) <C(\\F\\ m{v) + \\D*e\\ m[T) ) 
<C(||F|| ff i ( r) + ||e|| ff -i/2 (r) ) 

by the above mentioned continuity of D* . Now by the unique solvability of (58), 
we also have the stability estimate ||e|| H- x / 2 (r) — C\\F\\ H i/2( T y This, together 
with the continuous imbedding of H 1 ^) into i? 1 / 2 (r) shows that 

Nk(r) < C\\F\\ m{T) . (59) 

Wc will use this regularity result in the proof of the next theorem. 

Theorem 3. There is an H > and a < 5 < 1 such that whenever the coarse 
grid meshsize h\ is less than H , the error reduction operator E of Algorithm 1 
to the indefinite acoustic single layer equation satisfies 



||E||a < 6. (60) 
Here, S is independent of the refinement level J. 

Proof. This proof proceeds by verifying (47) and applying Theorem 2. To verify 
(47), we begin with the following: 

\V(a,fx)\ = \{Da,fi)\ < ||!><r|| Hl( r)||At|U-i(r) 

< C IMU-v^rollMllff-Hr)- ( 61 ) 

This is a consequence of Lemma 5. We shall use (61) several times below. 

Wc first prove (47) for k > 1. Define D k = V k P k - A k P k . Then (D k a, fi k ) = 
T>(a, Hk) for all a in i?~ 1 / 2 (r) and all fi k in M k . 



Z k = T k -T k = R k [V k P k - A k P k j = R k D k . 
For any a G Mj and k > 1, we have 

\\Z k a\\ 2 A = A{R k b k o,R k b k o) < \ k [R k D k a, R k D k a] k by (31) 
= A fc ^-{D k a, R k D k a)- X = V{a, R k D k a) by (29) 

< C ll cr llH-i/2(r)Pfc£>fecr||H-i(r) by (61). 

The last factor can be estimated by 

\\R k D k a\\%. 1{r) = (R k D h a,R k D h a)-! < C [R k D k o, R k D k o] k 

= C 2-(D k a,R k D k cT)-i< C ^\\D k o\\ H - Hr) \\R k D k a\\ H - Hr) , 
Afc \ k 
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and also noting that 

\\D k (r\\ 2 -i = (D k (r,Dko)-i = ([V h P k -A k P k )a,b h a 
= V{a,b k a) < C||a|| H -i/ 2(r) ||£) fc a|| H - 1( r). 
In combination, these show 

\\Zk<r\\A= \\RkD k cr\\A < C A^ 1/2 ||cr|| ff -i/2 (r) < C h X k /2 \\a\\ H -i / 2 ( Yy 
The last inequality follows from (31) and (38). Hence, for k > 2, 

\\Z k \\ = sup — - — < Ch k < , 

oeM h ||0l|A 

so we have verified (47) with a = 1/2. 

To prove (47) on the coarsest level (k = 1), we will use (59) and the following 
duality argument along the lines of a similar argument in [11]. Let a in AA k - 
Define 

This is a continuous linear functional on i/ _1 / 2 (r) and hence there is a unique 
solution e to (58) with this F. Hence, 



\\(I-P 1 )a\\H-i<T)=F((I-Pi)°) 
= V((I-P 1 )a,£) 
= V((I-P 1 )a,e-e 1 ) 

< C\\(I - Fi)cr|| H -i/ a(r) ||e - ei\\ H -^ r) (62) 

for any e\ in Mi. We choose an z\ with optimal approximation properties. 
Note that since 

\\F\\m { Y) = - PiMh-ht), 
the regularity result (59) holds for e. Therefore, 

||e - eiHa-i/^) < Ch\ /2 \\e\\ LHr) < Ch\ /2 \\F\\ Hl{r) 
= C^ /2 ||(/-P 1 ) < 7|| ff - 1(r) . 
Using this in (62), we conclude that 

||(7 - Pi)<7|| H - 1(r) < Ch\ /2 \\(I - PM\h-^ { t Y (63) 
We use (63) to estimate the norm of Z\ as follows. 

A(Z 1 a,n 1 )=A((P 1 -P 1 )a, m) = D {{I - Pj)a , Ml ) 

<a||(/-PxV||H-x ( r)||Mi||A by (61) 

< Ch\ /2 \\{I - Pi)<7|| H -i/ a(r) ||mi||a by (63) 

< C /i 1 //2 j|o'|| H -i/2( r -) ||/Lh||a by Lemma 3. 
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This proves (47) for the k = 1 case as well. 

Hence, we can apply Theorem 2 with a = 1/2 to get 

||E|| A < ||E|| A + C 2 /4 /2 < ~6 + C 2 h\ / \ 

where S is a positive number less than one given by Theorem 1. It is now clear 
that when hi is small enough, the result follows. 

4 Numerical Experiments 

4.1 Effect of weaker inner product on eigenfunctions 

The performance of the multigrid algorithm described in Section 2 depends 
crucially on the spectral behavior of the positive definite operator A. For a 
discrctized version of this integral operator the eigenfunctions corresponding to 
small magnitude eigenvalues are highly oscillatory, while those eigenfunctions 
corresponding to the large end of the spectrum are non-oscillatory. Standard 
multigrid approaches are successful for operator equations with the opposite 
spectral behavior and the use of the weaker inner products from Section 2.1 
effectively transforms the single layer problem into this setting. In this section 
we present the details of two examples showing the undesirable behaviour of the 
discrete eigenfunctions of the stiffness matrix associated with A while working in 
the (natural) H' 1 / 2 ^) inner product, and the effect of working with the weaker 
inner product instead. We demonstrate the effectiveness of this approach for 
two geometries, a (smooth) circle and a Lipschitz domain (square). 

Recall again that if the boundary Y is discrctized by means of a partition 
sci, %2 j ■ • ■ j xn, a^iv+i — x\, we denote by 4>i and U the characteristic function and 
respectively the length of the element Tj = conv(xj, Xi + \). The span of the {<fii} 
is denoted by Mk- In the case of the square, the boundary discretization consists 
of straight line segments; whereas in the case of the circle our discretization 
consists of arcs of equal angle. The basis coefficients of an element a € Mk 
with respect to the {4>i} are given in terms of the vector e(er) defined in (10). 

For a given discretization we are interested in the spectrum of the N x N 
stiffness matrix corresponding to the single layer operator for the Laplacian with 
entries [h\i,j = (A<j)j , <f>i). In Figure 1 the left-hand plots show the eigenfunction 
of A for a circular domain, corresponding to the smallest (top left) and largest 
(bottom left) eigenvalues respectively. These figures illustrate the phenomenon 
described above. 

In Section 2.1 we introduced a discrete inner product on the discretization 
space and proved that its associated norm is equivalent to the natural H~ 1 (T) 
norm. As before we denote by A the finite difference matrix corresponding to 
— u" + u on r with periodic boundary conditions, and by H the diagonal matrix 
with i th diagonal entry ij. According to the definition of the discrete inner 
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product in (12) we find for a, 6 £ A4k- 



[a,6] k = [e^-V)], , 0i>r = £ W^Mli ([«) ^ 

t j i 

= [e(^)rHA- 1 [e(a)] ! 

and also, 

A(0,0) = $>(0)] 4 PM A(&,&) = [e(0)]* A [e(0)]. 

The smoothing procedure, defined in Section 2.2, depends on the largest eigen- 
value of the following Rayleigh quotient with respect to the weaker inner prod- 
uct. 

A(0,0) e(0)*Ae(0) y* A*A A y 

A = SUP -: : = SUD — — : — — = SUP . 

eeM k [0,0} k eeM k e(0)* H A -1 e(0) y£ £ N y*k*Ey 

The two matrices A* A A and A* H are Hermitian. Therefore, A is the largest 
generalized eigenvalue of the problem: 

A* A A y = A A* H y, or cquivalently, A A y = A H y. (64) 

The corresponding cigenfunction is given in terms of its basis coefficients e(0) = 
A y. We note that the L 2 norm of is easy to compute. 

= E^M^W (^>=e(0)*He(0) = y*A*HAy. 

This allows us to normalize the coefficient vector e(0) such that has norm one 
in L 2 (r). We are now ready to examine the spectral behavior in the smoothing 
operation in terms of the generalized eigenvalue problem (64). 

In the case where T is a circle centered at the origin, the entries of the matrix 
A are particularly simple to compute. We require the radius to be bounded by 
R < 1/2 in order to guarantee the positive definiteness of the integral operator 
A and we discretize the circle with arcs Tj, i = 1..N of equal angle. 



[% = (A$,-,$i) = y J \og(\x - y\)ds x ds v 

■ J J log(2i? 2 (1 - cos(0, - Oj))) R 2 dOjdO i 



1 1 

'2tt 2 

47T 



r 4 Jt 



log(2R 2 ) dOjdOi 



R? 

47T 



log(l - cos(0 ?; - Oj)) dd ] d9 l 



hi 



47T 



log(2i? 2 ) 



2 R 2 e>2 f 



(65) 
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Hence, we have to evaluate the weakly singular integrals 

Iji = J J log (1 - cos(e t - Bj))) d6 3 d9 t 

for all possible choices of I\ and 1^-. Due to symmetry of the circle, the single 
layer matrix is a symmetric Toeplitz matrix and hence it is sufficient to compute 
its first row. For \i — j\ > 1, the integrand in Iij is smooth, and we use Gaussian 
quadrature to compute the entries. In the two remaining cases when either 
\i — j\ = 1 mod N or i = j, we separate the singularity of the integrand 
according to 

t 2 t 4 t 6 t 8 

log(l-COS(t))=log(-- - + --- + ...) 

= log(i 2 ) - log(2) + log(l - %- + % + . . .) (66) 



= fs(t) + fa(t). 

We then integrate the singular term f s (t) exactly and use Gaussian quadrature 
to compute the integral over the nonsingular function f a (t). 

Figure 1 shows the effect of the smoothing operator. The eigenfunction of 
the generalized eigenvalue problem AAy = XHy corresponding to the smallest 
eigenvalue (top right) is smooth (and represents a function of period 2t: on 
the circle. The eigenfunction of the largest eigenvalue (bottom right) is highly 
oscillatory. All eigenfunctions shown correspond to N = 300 elements on the 
circle. 

As a second test case we consider a square with side length 1/2. We compute 
the eigenfunctions corresponding to the four smallest and the four largest eigen- 
values for both the stiffness matrix A and the generalized eigenvalue problem 
AAy = XRy. Figures 2 and 3 illustrate the results for a uniform discretization 
of the boundary courve T with meshsize h = 1/50. Again, we observe that the 
eigenfunctions of the generalized eigenvalue problem display the reversed (and 
sought-after) smoothness behavior. The same behavior has been observed in 
cases of quasi-uniform meshes. 

4.2 Multigrid convergence results 

In this section we present numerical convergence results for the multigrid al- 
gorithm described in Section 2 to underline its effective use. Various results 
from a different set of experiments have already been reported in [13]. Recall 
that we want to solve the exterior Helmholtz problem with prescribed Dirich- 
let data on the boundary of a scattering object. The fundamental solution 
solves the Helmholtz equation away from its singularity x* 
and satisfies the correct growth condition at infinity. If we place the singularity 
into the interior of the scattering domain, we can use this so called point source 
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Figure 1: Eigcnfunction behavior, circular curve 
Eigcnfunctions for a circular curve with N = 300. They correspond to 
the smallest e-val of A (top left), the smallest c-val of hky = XEy (top right), 
the largest e-val of A (bottom left) and the largest c-val of AAy = XUy (bottom 
right). 




Figure 2: Eigcnfunction behavior, square boundary curve, small eigenvalues 
Eigenfunction behavior for a square boundary curve T and a uniform mesh with 
mesh size h = 1/50. Eigenfunctions correspond to the four smallest 
e-vals of A (left) and the four smallest generalized e-vals of hky = \Ry (right). 
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Figure 3: Eigenfunction behavior, square boundary curve, large eigenvalues 
Eigenfunction behavior for a square boundary curve T and a uniform mesh with 
mesh size h = 1 /50. Eigenfunctions correspond to the four largest 
e-vals of A (left) and the four largest generalized e-vals of Aky = AHy (right). 

problem as a convenient test case, for which the exact solution is known. We 
also present results for the scattering of an incident plane wave by polygonal 
shaped obstacles. In all the tables below, H is the coarsest mesh size, and h is 
the finest mesh. 

The fast implementation of the matrix-vector multiplications in the algo- 
rithm was not subject of our study and hence we do not report on the over- 
all CPU time used by the algorithm. We anticipate that the CPU time will 
be competitive once those matrix operations are implemented using matrix- 
compression techniques such as H-matrices. 

4.2.1 Effect of domain shape on performance 

We first present results for point-source scattering from 4 different objects: a 
square, a rectangle (sides of ratio 1:4), an equilateral triangle and a thin wedge. 
The thin wedge is described in terms of the xy-coordinatcs of its three corner 
points, namely (0,0), (l/3,v / 3/3) and (0,1/15). This amount to an angle of 
7r/3 between the horizontal x-axis and the lower edge of the wedge. In each 
of these examples, the diameter of the object is less than 1, which guarantees 
the positive definiteness of the potential single layer operator in the sense of 
(4). In each of Tables 1-4, the proposed multigrid scheme is used as a linear 
solver. We report the number of iteration numbers required to reach a given 
relative residual norm. The point source is located inside the domain, so the 
true solution is known in each case. 
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Table 1: Linear multigrid iteration counts with k — 2.1, point-source inside 
square, 10 ~ 6 relative residual norm 





1/2 


H 

1/4 


1/8 


1/16 


1/32 


Degrees 
of freedom 


1/4 


17 










32 


1/8 


16 


15 








64 


1/16 


15 


15 


15 






128 


h 1/32 


15 


15 


15 


15 




256 


1/64 


16 


16 


15 


15 


15 


512 


1/128 


16 


16 


16 


16 


16 


1024 


1/256 


16 


16 


16 


16 


16 


2048 



Table 2: Linear multigrid iteration counts with K = 2.1, point-source inside 
rectangle, 10 -6 relative residual norm 





1/2 


H 

1/4 


1/8 


1/16 


1/32 


Degrees 
of freedom 


1/4 


15 










72 


1/8 


17 


15 








144 


1/16 


18 


18 


16 






288 


h 1/32 


18 


19 


18 


16 




576 


1/64 


19 


19 


19 


18 


16 


1152 


1/128 


19 


19 


19 


19 


19 


2304 



Table 3: Linear multigrid iteration counts with n = 2.1, point-source inside 
triangle, 10~ 6 relative residual norm 





1 


1/2 


1/4 


H 

1/8 


1/16 


1/32 


1/64 


1/128 


Dcgr. of 
freedom 


1/2 


21 
















6 


1/4 


33 


30 














12 


1/8 


37 


37 


36 












24 


1/16 


28 


28 


27 


27 










48 


h 1/32 


24 


24 


24 


23 


23 








96 


1/64 


23 


23 


23 


22 


22 


22 






192 


1/128 


22 


22 


22 


22 


22 


21 


22 




384 


1/256 


21 


21 


21 


21 


21 


21 


21 


21 


768 


1/512 


20 


20 


20 


21 


21 


21 


20 


20 


1536 
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Table 4: Linear multigrid iteration counts with k = 27r, thin wedge, pre- 
smoothing only, 10 -4 relative residual norm 





1/4 


H 

1/8 


1/16 


1/32 


1/64 


Degr. of 
freedom 


1/16 


31 


21 








352 


h 1/32 


18 


13 


9 






704 


1/64 


17 


11 


8 


5 




1408 


1/128 


17 


11 


7 


4 


3 


2816 



Table 5: Linear multigrid iteration counts with k = 10.2, point source inside 
triangle, pre-smoothing only, 10 -6 relative residual norm 





1/4 


1/8 


H 

1/16 


1/32 


1/64 


1/128 


Degr. of 
freedom 


1/8 


38 












24 


1/16 


36 


33 










48 


h 1/32 


32 


30 


28 








96 


1/64 


27 


26 


25 


24 






192 


1/128 


24 


23 


22 


22 


22 




384 


1/256 


22 


21 


21 


21 


21 


22 


768 


1/512 


21 


21 


21 


21 


21 


21 


1536 



4.2.2 Multigrid as linear solver/ preconditioner 

In this section, we provide results about the use of the proposed multigrid 
scheme as a linear solver (Tables ??), and as a preconditioner for GMRES 
(used without restart, Tables ??). We present iteration counts in each case 
with just presmoothing, or with both pre-and post-smoothing. We use the 
same equilateral triangle as the domain as in the previous section. Again, a 
point-source is placed inside the domain, so we can compare with the exact 
solution. In contrast to the previous section, here we present results for wave 
number k = 10.2. The tolerances in relative residual norm are 10~ 6 , 10~ 9 for the 
linear solver and preconditioned GMRES, respectively. For both, multigrid as 
a linear solver and preconditioned GMRES, the iteration numbers stay almost 
constant with increasing degress of freedom, whereas the iteration numbers of 
GMRES grow considerably. Convergence for the algorithm with an additional 
postsmoothing step follows from our convergence result by standard arguments 
in multigrid theory. 

4.2.3 Effect of frequency on performance 

Here we present the effect of increasing k on the number of multigrid iterations 
taken to achieve a given relative residual error. We place a point source inside 
a square whose diameter is less than 1. We expect the coarsest mesh required 
should satisfy the constraint kH w constant. At least in this example, the 
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Table 6: Linear multigrid iteration counts with k = 10. 2, point source inside 
triangle, pre- and post-smoothing, 1CP 6 relative residual norm 





1/4 


1/8 


H 

1/16 


1/32 


1/64 


1/128 


Dcgr. of 
freedom 


1/8 


19 












24 


1/16 


20 


17 










48 


h 1/32 


18 


16 


14 








96 


1/64 


17 


15 


14 


13 






192 


1/128 


15 


14 


13 


13 


13 




384 


1/256 


15 


13 


13 


13 


13 


13 


768 


1/512 


16 


14 


13 


13 


13 


13 


1536 



Table 7: GMRES iteration counts, K = 10. 2, point source inside triangle, trian- 



gle, pre-smoothing only, 10 9 relative residual norm 





H 






GMRES 




1/2 


1/4 


1/8 


1/16 


without preconditioning 


1/8 


17 


16 






16 


1/16 


23 


21 


20 




24 


h 1/32 


23 


22 


22 


19 


31 


1/64 


24 


23 


23 


22 


37 


1/128 


25 


25 


24 


23 


44 


1/256 


25 


27 


24 


23 


52 


1/512 


26 


28 


24 


24 


63 


1/1024 


27 


28 


26 


26 


74 



Table 8: GMRES iteration counts, K = 10.2, triangle, uniform grid, 6 = 1, pre 



and post-smoothing, 10 9 relative residual norm 





H 






GMRES 




1/2 


1/4 


1/8 


1/16 


without preconditioning 


1/8 


18 


15 






16 


1/16 


20 


17 


17 




242 


h 1/32 


20 


18 


18 


17 


31 


1/64 


21 


19 


19 


18 


37 


1/128 


21 


21 


19 


19 


44 


1/256 


22 


22 


21 


20 


52 


1/512 


22 


24 


21 


21 


63 


1/1024 


23 


24 


21 


21 


74 
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Table 9: Linear multigrid iteration counts with k = 2.1, point source inside 
square, pre- and post-smoothing, 10~ 6 relative residual norm 



H 





1/4 


1/8 


1/16 


1/32 


1/64 


1/128 


1/32 


12 


12 


12 








1/64 


13 


13 


12 


12 






1/128 


13 


13 


13 


13 


12 




1/256 


13 


13 


13 


13 


13 


13 


1/512 


14 


14 


14 


14 


14 


13 



1/256 



13 



Degr. of 
freedom 
256 
512 
1024 
2048 
4096 



Table 10: Linear multigrid iteration counts with k = 10. 2, point source inside 



square, pre- and post-smoothing, 10 6 
not converge for H=l/4. 

H 



relative residual norm. The method did 





1/4 


1/8 


1/16 


1/32 


1/64 


1/128 


1/32 


* 


15 


13 








1/64 


* 


14 


12 


11 






1/128 


* 


14 


12 


11 


11 




1/256 


* 


13 


12 


11 


11 


11 


1/512 


* 


12 


12 


12 


12 


12 



1/256 



12 



Degr. of 
freedom 
256 
512 
1024 
2048 
4096 



method performs well even though this constraint was not strictly satisfied. For 
example, in Table 9, kH = 0.252, while in Table 11, kH = 1.575, 

The numerical experiments show that the multigrid algorithm presented 
and analyzed in this paper is an efficient tool to solve the first kind single layer 
equation when used as a linear solver or as a preconditioning procedure for other 
solvers such as GMRES. 



Table 11: Linear multigrid iteration counts with n = 50. 4, point source inside 



square, pre- and post-smoothing, 10 6 
not converge for H> 1/32. 

H 



relative residual norm. The method did 





1/4 


1/8 


1/16 


1/32 

* 


1/64 


1/128 


1/32 


* 


* 






1/64 


* 


* 


* 


* 






1/128 


* 


* 




* 


13 




1/256 


* 


* 




* 


13 


11 


1/512 


* 


* 




* 


12 


11 



1/256 



11 



Degr. of 
freedom 
256 
512 
1024 
2048 
4096 
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Figure 4: Plane-wave scattering from wedge: linear MG is used to compute the 
unknown density, and numerical quadrature is used to reconstruct the scattered 
field. 
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